!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine eigenspace (eigenvalue, eigenvector_l, eigenvector_r, ro, mx, my, mz, bx, by, bz, en, nx, gm)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx
   real(8), intent(in) :: gm
   real(8), dimension(nx), intent(in) :: ro, mx, my, mz, bx, by, bz, en
   real(8), dimension(8, nx), intent(inout) :: eigenvalue
   real(8), dimension(8, 8, nx), intent(inout) :: eigenvector_l, eigenvector_r

   integer(4) :: i
   real(8) :: vx, vy, vz, b2, v2, pr, c, c2, c22, s2, cax, ca2, cfx, csx, cfx2, csx2, eps, gmm1, ro_sqrt
   real(8) :: sgn, beta_y, beta_z, alpha_f, alpha_s, gm_1, gm_2, tao, gm_f, gm_a, gm_s

!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12
   gmm1 = gm - 1.0d0

   do i = 1, nx
      vx = mx(i) / (ro(i) + eps)
      vy = my(i) / (ro(i) + eps)
      vz = mz(i) / (ro(i) + eps)
      b2 = bx(i) * bx(i) + by(i) * by(i) + bz(i) * bz(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i) - v2 * ro(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr
      s2 = c2 + b2
      ca2 = bx(i) * bx(i)
      cfx2 = (s2 + sqrt (s2 * s2 - 4.0d0 * c2 * ca2)) / 2.0d0
      csx2 = (s2 - sqrt (s2 * s2 - 4.0d0 * c2 * ca2)) / 2.0d0
      if (isnan(cfx2) .or. cfx2 .gt. 1.0d0 / eps .or. cfx2 .lt. -1.0d0 / eps) cfx2 = eps
      if (isnan(csx2) .or. csx2 .gt. 1.0d0 / eps .or. csx2 .lt. -1.0d0 / eps) csx2 = eps

      if (bx(i) .ge. 0.0d0) then
         sgn = 1.0d0
      else
         sgn = -1.0d0
      endif

      b2 = by(i) * by(i) + bz(i) * bz(i)

      if (b2 .gt. eps) then
         beta_y = by(i) / sqrt(b2)
         beta_z = bz(i) / sqrt(b2)
      else
         beta_y = 1.0d0 / sqrt(2.0d0)
         beta_z = 1.0d0 / sqrt(2.0d0)
      endif

      if (b2 > eps .or. abs (gm * pr - bx(i) * bx(i) .gt. eps)) then
         alpha_f = sqrt (c2 - csx2) / sqrt (cfx2 - csx2)
         alpha_s = sqrt (cfx2 - c2) / sqrt (cfx2 - csx2)
      else
         alpha_f = 1.0d0 / sqrt(2.0d0)
         alpha_s = 1.0d0 / sqrt(2.0d0)
      endif

      c = sqrt (c2 / ro(i))
      cax = sqrt (ca2 / ro(i))
      cfx = sqrt (cfx2 / ro(i))
      csx = sqrt (csx2 / ro(i))
      ro_sqrt = sqrt (ro(i))
      c22 = c * c * 2.0d0

      eigenvalue(1, i) = vx - cfx
      eigenvalue(2, i) = vx - cax
      eigenvalue(3, i) = vx - csx
      eigenvalue(4, i) = vx
      eigenvalue(5, i) = vx
      eigenvalue(6, i) = vx + csx
      eigenvalue(7, i) = vx + cax
      eigenvalue(8, i) = vx + cfx

      gm_1 = gmm1 / 2.0d0
      gm_2 = (gm - 2.0d0) / gmm1
      tao = gmm1 / c2 * ro(i)
      gm_f = alpha_f * cfx * vx - alpha_s * csx * sgn * (beta_y * vy + beta_z * vz)
      gm_a = sgn * (beta_z * vy - beta_y * vz)
      gm_s = alpha_s * csx * vx + alpha_f * cfx * sgn * (beta_y * vy + beta_z * vz)

      eigenvector_r(1, 1, i) = alpha_f
      eigenvector_r(1, 2, i) = alpha_f * (vx - cfx)
      eigenvector_r(1, 3, i) = alpha_f * vy + csx * alpha_s * beta_y * sgn
      eigenvector_r(1, 4, i) = alpha_f * vz + csx * alpha_s * beta_z * sgn
      eigenvector_r(1, 5, i) = 0.0d0
      eigenvector_r(1, 6, i) = c * alpha_s * beta_y / ro_sqrt
      eigenvector_r(1, 7, i) = c * alpha_s * beta_z / ro_sqrt
      eigenvector_r(1, 8, i) = alpha_f * (v2 / 2.0d0 + cfx * cfx - gm_2 * c * c) - gm_f
      eigenvector_r(8, 1, i) = alpha_f
      eigenvector_r(8, 2, i) = alpha_f * (vx + cfx)
      eigenvector_r(8, 3, i) = alpha_f * vy - csx * alpha_s * beta_y * sgn
      eigenvector_r(8, 4, i) = alpha_f * vz - csx * alpha_s * beta_z * sgn
      eigenvector_r(8, 5, i) = 0.0d0
      eigenvector_r(8, 6, i) = c * alpha_s * beta_y / ro_sqrt
      eigenvector_r(8, 7, i) = c * alpha_s * beta_z / ro_sqrt
      eigenvector_r(8, 8, i) = alpha_f * (v2 / 2.0d0 + cfx * cfx - gm_2 * c * c) + gm_f
      eigenvector_r(2, 1, i) = 0.0d0
      eigenvector_r(2, 2, i) = 0.0d0
      eigenvector_r(2, 3, i) = -beta_z * sgn
      eigenvector_r(2, 4, i) = beta_y * sgn
      eigenvector_r(2, 5, i) = 0.0d0
      eigenvector_r(2, 6, i) = -beta_z / ro_sqrt
      eigenvector_r(2, 7, i) = beta_y / ro_sqrt
      eigenvector_r(2, 8, i) = -gm_a
      eigenvector_r(7, 1, i) = 0.0d0
      eigenvector_r(7, 2, i) = 0.0d0
      eigenvector_r(7, 3, i) = -beta_z * sgn
      eigenvector_r(7, 4, i) = beta_y * sgn
      eigenvector_r(7, 5, i) = 0.0d0
      eigenvector_r(7, 6, i) = beta_z / ro_sqrt
      eigenvector_r(7, 7, i) = -beta_y / ro_sqrt
      eigenvector_r(7, 8, i) = -gm_a
      eigenvector_r(3, 1, i) =  alpha_s
      eigenvector_r(3, 2, i) =  alpha_s * (vx - csx)
      eigenvector_r(3, 3, i) =  alpha_s * vy - cfx * alpha_f * beta_y * sgn
      eigenvector_r(3, 4, i) =  alpha_s * vz - cfx * alpha_f * beta_z * sgn
      eigenvector_r(3, 5, i) =  0.0d0
      eigenvector_r(3, 6, i) =  -c * alpha_f * beta_y / ro_sqrt
      eigenvector_r(3, 7, i) =  -c * alpha_f * beta_z / ro_sqrt
      eigenvector_r(3, 8, i) =  alpha_s * (v2 / 2.0d0 + csx * csx - gm_2 * c * c) - gm_s
      eigenvector_r(6, 1, i) =  alpha_s
      eigenvector_r(6, 2, i) =  alpha_s * (vx + csx)
      eigenvector_r(6, 3, i) =  alpha_s * vy + cfx * alpha_f * beta_y * sgn
      eigenvector_r(6, 4, i) =  alpha_s * vz + cfx * alpha_f * beta_z * sgn
      eigenvector_r(6, 5, i) =  0.0d0
      eigenvector_r(6, 6, i) =  -c * alpha_f * beta_y / ro_sqrt
      eigenvector_r(6, 7, i) =  -c * alpha_f * beta_z / ro_sqrt
      eigenvector_r(6, 8, i) =  alpha_s * (v2 / 2.0d0 + csx * csx - gm_2 * c * c) + gm_s
      eigenvector_r(4, 1, i) =  1.0d0
      eigenvector_r(4, 2, i) =  vx
      eigenvector_r(4, 3, i) =  vy
      eigenvector_r(4, 4, i) =  vz
      eigenvector_r(4, 5, i) =  0.0d0
      eigenvector_r(4, 6, i) =  0.0d0
      eigenvector_r(4, 7, i) =  0.0d0
      eigenvector_r(4, 8, i) =  v2 / 2.0d0
      eigenvector_r(5, 1, i) =  0.0d0
      eigenvector_r(5, 2, i) =  0.0d0
      eigenvector_r(5, 3, i) =  0.0d0
      eigenvector_r(5, 4, i) =  0.0d0
      eigenvector_r(5, 5, i) =  1.0d0
      eigenvector_r(5, 6, i) =  0.0d0
      eigenvector_r(5, 7, i) =  0.0d0
      eigenvector_r(5, 8, i) =  bx(i)

      eigenvector_l(1, 1, i) =  (gm_1 * alpha_f * v2 + gm_f) / c22
      eigenvector_l(2, 1, i) =  (-gmm1 * alpha_f * vx - alpha_f * cfx) / c22
      eigenvector_l(3, 1, i) =  (-gmm1 * alpha_f * vy + csx * alpha_s * beta_y * sgn) / c22
      eigenvector_l(4, 1, i) =  (-gmm1 * alpha_f * vz + csx * alpha_s * beta_z * sgn) / c22
      eigenvector_l(5, 1, i) =  (-gmm1 * alpha_f * bx(i)) / c22
      eigenvector_l(6, 1, i) =  (-gmm1 * alpha_f * by(i) + ro_sqrt * c * alpha_s * beta_y) / c22
      eigenvector_l(7, 1, i) =  (-gmm1 * alpha_f * bz(i) + ro_sqrt * c * alpha_s * beta_z) / c22
      eigenvector_l(8, 1, i) =  gmm1 * alpha_f / c22
      eigenvector_l(1, 8, i) =  (gm_1 * alpha_f * v2 - gm_f) / c22
      eigenvector_l(2, 8, i) =  (-gmm1 * alpha_f * vx + alpha_f * cfx)/c22
      eigenvector_l(3, 8, i) =  (-gmm1 * alpha_f * vy - csx * alpha_s * beta_y * sgn) / c22
      eigenvector_l(4, 8, i) =  (-gmm1 * alpha_f * vz - csx * alpha_s * beta_z * sgn) / c22
      eigenvector_l(5, 8, i) =  (-gmm1 * alpha_f * bx(i)) / c22
      eigenvector_l(6, 8, i) =  (-gmm1 * alpha_f * by(i) + ro_sqrt * c * alpha_s * beta_y) / c22
      eigenvector_l(7, 8, i) =  (-gmm1 * alpha_f * bz(i) + ro_sqrt * c * alpha_s * beta_z) / c22
      eigenvector_l(8, 8, i) =  gmm1 * alpha_f / c22
      eigenvector_l(1, 2, i) =  gm_a / 2.0d0
      eigenvector_l(2, 2, i) =  0.0d0
      eigenvector_l(3, 2, i) =  -beta_z * sgn / 2.0d0
      eigenvector_l(4, 2, i) =  beta_y * sgn / 2.0d0
      eigenvector_l(5, 2, i) =  0.0d0
      eigenvector_l(6, 2, i) =  -ro_sqrt * beta_z / 2.0d0
      eigenvector_l(7, 2, i) =  ro_sqrt * beta_y / 2.0d0
      eigenvector_l(8, 2, i) =  0.0d0
      eigenvector_l(1, 7, i) =  gm_a / 2.0d0
      eigenvector_l(2, 7, i) =  0.0d0
      eigenvector_l(3, 7, i) =  -beta_z * sgn / 2.0d0
      eigenvector_l(4, 7, i) =  beta_y * sgn / 2.0d0
      eigenvector_l(5, 7, i) =  0.0d0
      eigenvector_l(6, 7, i) =  ro_sqrt * beta_z / 2.0d0
      eigenvector_l(7, 7, i) =  -ro_sqrt * beta_y / 2.0d0
      eigenvector_l(8, 7, i) =  0.0d0
      eigenvector_l(1, 3, i) =  (gm_1 * alpha_s * v2 + gm_s) / c22
      eigenvector_l(2, 3, i) =  (-gmm1 * alpha_s * vx - alpha_s * csx) / c22
      eigenvector_l(3, 3, i) =  (-gmm1 * alpha_s * vy - cfx * alpha_f * beta_y * sgn) / c22
      eigenvector_l(4, 3, i) =  (-gmm1 * alpha_s * vz - cfx * alpha_f * beta_z * sgn) / c22
      eigenvector_l(5, 3, i) =  (-gmm1 * alpha_s * bx(i)) / c22
      eigenvector_l(6, 3, i) =  (-gmm1 * alpha_s * by(i) - ro_sqrt * c * alpha_f * beta_y) / c22
      eigenvector_l(7, 3, i) =  (-gmm1 * alpha_s * bz(i) - ro_sqrt * c * alpha_f * beta_z) / c22
      eigenvector_l(8, 3, i) =  gmm1 * alpha_s / c22
      eigenvector_l(1, 6, i) =  (gm_1 * alpha_s * v2 - gm_s) / c22
      eigenvector_l(2, 6, i) =  (-gmm1 * alpha_s * vx + alpha_s * csx) / c22
      eigenvector_l(3, 6, i) =  (-gmm1 * alpha_s * vy + cfx * alpha_f * beta_y * sgn) / c22
      eigenvector_l(4, 6, i) =  (-gmm1 * alpha_s * vz + cfx * alpha_f * beta_z * sgn) / c22
      eigenvector_l(5, 6, i) =  (-gmm1 * alpha_s * bx(i)) / c22
      eigenvector_l(6, 6, i) =  (-gmm1 * alpha_s * by(i) - ro_sqrt * c * alpha_f * beta_y) / c22
      eigenvector_l(7, 6, i) =  (-gmm1 * alpha_s * bz(i) - ro_sqrt * c * alpha_f * beta_z) / c22
      eigenvector_l(8, 6, i) =  gmm1 * alpha_s / c22
      eigenvector_l(1, 4, i) =  1.0d0 - tao * v2 / 2.0d0
      eigenvector_l(2, 4, i) =  vx * tao
      eigenvector_l(3, 4, i) =  vy * tao
      eigenvector_l(4, 4, i) =  vz * tao
      eigenvector_l(5, 4, i) =  tao * bx(i)
      eigenvector_l(6, 4, i) =  by(i) * tao
      eigenvector_l(7, 4, i) =  bz(i) * tao
      eigenvector_l(8, 4, i) =  -tao
      eigenvector_l(1, 5, i) =  0.0d0
      eigenvector_l(2, 5, i) =  0.0d0
      eigenvector_l(3, 5, i) =  0.0d0
      eigenvector_l(4, 5, i) =  0.0d0
      eigenvector_l(5, 5, i) =  1.0d0
      eigenvector_l(6, 5, i) =  0.0d0
      eigenvector_l(7, 5, i) =  0.0d0
      eigenvector_l(8, 5, i) =  0.0d0

   enddo

   return
end subroutine eigenspace
